Zhao Liang, Xu Shun, Tu Yu-Song, Zhou Xin. Determination of the vapor–liquid transition of square-well particles using a novel generalized-canonical-ensemble-based method. Chinese Physics B, 2017, 26(6): 060202
Permissions
Determination of the vapor–liquid transition of square-well particles using a novel generalized-canonical-ensemble-based method
Zhao Liang1, Xu Shun2, Tu Yu-Song1, †, Zhou Xin3, ‡
College of Physical Science and Technology, Yangzhou University, Yangzhou 225009, China
Supercomputer Center, Computer Network Information Center, Chinese Academy of Sciences, Beijing 100190, China
School of Physical Sciences, University of Chinese Academy of Sciences, Beijing 100049, China
Project supported by the National Natural Science Foundation for Outstanding Young Scholars, China (Grant No. 11422542), the National Natural Science Foundation of China (Grant Nos. 11605151 and 11675138), and the Shanghai Supercomputer Center of China and Special Program for Applied Research on Super Computation of the NSFC–Guangdong Joint Fund (the second phase).
Abstract
The square-well (SW) potential is one of the simplest pair potential models and its phase behavior has been clearly revealed, therefore it has become a benchmark for checking new theories or numerical methods. We introduce the generalized canonical ensemble (GCE) into the isobaric replica exchange Monte Carlo (REMC) algorithm to form a novel isobaric GCE-REMC method, and apply it to the study of vapor–liquid transition of SW particles. It is validated that this method can reproduce the vapor–liquid diagram of SW particles by comparing the estimated vapor–liquid binodals and the critical point with those from the literature. The notable advantage of this method is that the unstable vapor–liquid coexisting states, which cannot be detected using conventional sampling techniques, are accessed with a high sampling efficiency. Besides, the isobaric GCE-REMC method can visit all the possible states, including stable, metastable or unstable states during the phase transition over a wide pressure range, providing an effective pathway to understand complex phase transitions during the nucleation or crystallization process in physical or biological systems.
A square-well (SW) potential is one of the simplest interaction models[1–4] and has attracted a huge variety of studies on its structural and thermodynamical properties during the past several decades.[5–15] The SW potential only consists of a hard sphere repulsion and a finite-range constant attraction, and the interaction between a pair of particles separated by a distance r is in the following form:
where σ is the hard-core diameter of the particle, λ is the range of the attraction potential, and is the attraction well depth. Despite its simplicity, the SW exhibits the vapor–liquid, vapor–solid, and liquid–solid phase transitions,[10] capturing the main phenomenology of real atomic fluids[16,17] even the biological solutions[18] and providing an insight into the understanding of complex physical and biological phenomena,[19–22] such as colloid nucleation and protein crystallization. To date, phase behavior of SW particles has been clearly revealed, therefore this model has turned into an important benchmark for checking new theories or numerical methods.[7,10,23–25]
Many theoretical or numerical methods have been applied to the exploration of the phase diagram of SW particles, such as Gibbs ensemble Monte Carlo (GEMC),[26] grand canonical Monte Carlo,[8] Gibbs-Duhem integration,[27] histogram reweighting,[28] discontinuous molecular dynamics,[7] and perturbation theory.[23] The most frequently used method is GEMC as coexistent binodals can be directly determined by setting two physically separated but thermodynamically equilibrated distinct coexistent bulk phases. However, the thermodynamical quantities of the critical point need to be extrapolated, and the high density phase is not efficiently studied because of the poor probability of particle insertion. Another popular method is the replica exchange Monte Carlo (REMC),[25,29] which involves the simultaneously parallel sampling in replicas of the same system at different temperatures and the exchange of configurations between replicas during the simulation process. It has been found to be effective in enhancing the sampling of the low-energy states on a free-energy landscape with many local minima, and thermodynamical information at different pressures or temperatures can be collected after running a single parallel simulation. Although this method accesses the phase equilibria of many single or multicomponent systems, configurations in the phase-coexisting regions will not be efficiently exchanged with each other as the sampling of configurations is hampered by the high energy barriers between the coexisting distinct phases. Despite these efforts in developing various methods to enhance the sampling of possible states during phase transitions, the difficulty in detecting the metastable states or unstable phase-coexisting states remains.
Recently, several novel algorithms beyond the traditional MC have been proposed, aiming to enhance the sampling efficiency or to predict the critical properties more precisely. A series of algorithms to enhance the sampling[30–33] of possible states has been systematically formulated based on a new ensemble, generalized canonical ensemble (GCE),[34–36] which overcomes the inability of sampling the metastable states or unstable phase-coexisting states during phase transitions by using the conventional numerical methods. They introduced the GCE into the canonical parallel tempering Monte Carlo (PTMC)[36] or the isobaric replica exchange molecular dynamics (REMD).[35] These simulation techniques with the combination of GCE, yield the unimodal Gaussian-like energy distributions with significant overlaps, even if when system states are in phase-coexisting regions, suggesting that by using these GCE-based methods, all the possible states even the unstable phase-coexisting states can be visited. Alternative classes of algorithms are designed beyond the detailed balance condition and also show the promising future for traditional MC.[37–39] The event-chain Monte Carlo,[40,41] has been proposed and applied to the simulation of hard-sphere with repulsive power-law interaction. They identified the first order liquid-hexatic and hexatic-solid transition. Similar irreversible tricks, designed by Hu et al., have been successfully generalized to the lattice system.[42] The critical point of self-avoiding walk was determined with approximately eight times more accuracy than available predictions in the literature.
In this work, we introduce the GCE into the REMC to form a novel GCE-based method, namely the isobaric GCE-REMC, then check its validity and further show advantages in determining the phase diagram of vapor–liquid transition of SW particles. The rest of this paper is organized as follows. Firstly, we briefly introduce the principle of GCE and derive the formula of acceptance probability for three types of attempting movements in the GCE-based isobaric REMC method. Secondly, this novel method is applied to the study of the vapor–liquid phase diagram of SW particles. The validity of this method is checked by calculating the vapor–liquid binodals and critical point, and also by comparing them with those from the literature. Finally, advantages in detecting all possible states during the vapor–liquid transition are demonstrated by analyzing the sampling efficiency and swap acceptance ratio between neighbor replicas.
2. Method and theory
2.1. Generalized canonical ensemble (GCE)
Simulation ensembles[43–45] are distinguished by selecting the effective potential U(r) in the conformation distribution function W(r),
where r denotes the atomic position vector of a conformation in the high-dimensional space. The corresponding energy distribution is
where S(E) is the micro-canonical entropy of the system (logarithm of the density of states in the energy E).
where β0, E0, and α are the referential inverse temperature, referential internal energy and control parameter, respectively. The energy distribution P(E) around the average energy in GCE is the approximate expression of the Gaussian-like form
indicating that the unstable coexistence and spinodal states with in the canonical ensemble () can thus be accessed in GCE by selecting .
2.2. Isobaric GCE-REMC algorithm
The isobaric GCE-REMC algorithm is similar to that of the normal isobaric REMC except replacing the internal energy by the GCE effective potential, which comprises three types of trial moves:
Displacement of a randomly selected particle.
In isobaric simulations, the internal energy E is replaced by the enthalpy H and the effective potential is with H = E + PV. Suppose that we have generated a trial move from state old to state new, the acceptance probability is
where the effective energy difference is
with the enthalpy difference .
Rescale of system volume.
To maintain the pressure at a fixed value, the system volume will be adjusted. In conventional isobaric simulations, the acceptance probability of volume changing from in state old to in state new is given by[46]
with the effective energy difference
In our isobaric GCE simulations, as the effective potential includes the square term , the effective energy difference is thus similar to the expression in Eq. (7) and is in the following form:
Exchange of system configuration and volume between neighbor replicas.
In isobaric GCE-REMC simulations, M replicas evolve from the same initial configuration in parallel. They have the same reference inverse temperatures and control parameters but different reference enthalpies. For convenience, we take , and evenly distributed Hi0 in the interval , . Each replica only exchanges its conformation with two neighbors, i.e., replicas can be classified into two sets: [(1, 2), (3, 4), …] or [(2, 3), (4, 5), …] and these two sets are randomly used. The number of replica pairs is . Supposing that inverse temperatures of replicas i, j are , , we have the acceptance probability of replica exchange
2.3. Determinations of binodals and the critical point
The visiting probability of an NPT system in the enthalpy space is
where β0 is the inverse value of . With the expression[36], the terms on the exponent can be further expressed as
In a first order transition system, the probability of finding the vapor state or the liquid state at the coexistent temperature T0 should be equal. This means that the area under the line or the curve , above the H axis ranging from Hl to Hv should be equal in the β–H plane as shown in Fig. 1. According to this equal area rule, the inverse temperature β0, enthalpy of vapor Hv and that of Hl in the coexistent vapor–liquid states can be determined. In addition, the spinodals are given under the condition of , i.e., the E and C points shown in Fig. 1.
Fig. 1. Schematic representation of determination of the vapor and liquid binodals. The inverse temperature β0 is determined by satisfying the equal area rule , i.e., the shadow areas above and under the line should be equal to each other, i.e., .
The critical point is recognized when the vapor and liquid binodals converge to a single point, where the enthalpies of vapor Hv and liquid Hl are equal to each other.
2.4. Simulation details
Isobaric GCE-REMC simulations were performed with M = 72 replicas at constant pressure. All systems evolve from the same initial configuration corresponding to a liquid state, where 300 particles are randomly dispersed in a cubic box with a dimension of 7.3 nm × 7.3 nm × 7.3 nm. Periodic boundary conditions are applied in all directions. The maximum particle displacement and volume rescale are adjusted to give ∼50% of acceptance for trial movements. The values of reference inverse temperature β0 of 72 replicas are 1.0 but the values of reference enthalpy H0 are taken from an evenly spaced interval [−1550, 50], and their values of sampling slope α are 0.006. There are 71 replica pairs and replica exchanges between neighbor replicas are performed every 100 MC cycles. The total MC cycles are 3 × 108 and data during the last 1 × 108 cycles are collected for further analyses.
3. Results and discussion
In the following discussion, we focus on the vapor–liquid transition of SW particles with the attraction range parameter λ = 1.5 at various pressures. It has been known that there exists an evident vapor–liquid coexistence line for ,[10] which will be convenient for checking the validity and show advantages of the novel isobaric GCE-REMC method. The reduced units , , , and will be used throughout the discussion.
Figure 2(a) shows the phase diagram of SW particles with the attraction range parameter λ = 1.5 in the density-temperature plane at various pressures. At pressure , isobaric lines all show similar oscillation behaviors with a relative minimum and a relative maximum, recognized as the vapor and liquid spinodals, respectively. The system state approaches to the vapor or liquid spinodals at constant pressure when the temperature of supercooled vapor continues to fall or that of superheated liquid continues to rise. Between these two spinodals, of the system state has a positive value, indicating that heating up the system will reduce the volume and increase the number density of particles, which is impossible to occur and violates the thermodynamical stability conditions. These states captured by our isobaric GCE-REMC simulations are in fact unstable, representing the phase-coexisting regions, and cannot be detected using the conventional simulation methods. As the pressure increases to , the oscillation disappears and the isobaric line shows a relatively flat plateau. The two vapor and liquid spinodals converge to a critical point with the estimated (, , , which compares favorably with the values of critical points from the literature, as listed in Table 1.
Fig. 2. (color online) (a) Phase diagram of SW particles with attraction range parameter λ = 1.5 in the temperature-density plane. (b) Comparisons of vapor–liquid binodals, critical points from our calculations and the literature. Binodals calculated at pressures of 0.08, 0.06, 0.04, and 0.02 are represented by solid symbols colored in orange, green, purple and blue, respectively. The crosses are taken from Vega et al.,[6] the pluses from Singh et al.,[9] and the empty squares from del Rio et al.[10] The vapor–liquid critical point in this work is estimated (, , , represented by the yellow square, in comparison with data from the literature represented by the black circles (for the details see Table 1).
Table 1.
Table 1.
Table 1.
Critical vapor–liquid point data of SW particles with λ = 1.5 from our isobaric GCE-REMC simulations and those from the literature.
Critical vapor–liquid point data of SW particles with λ = 1.5 from our isobaric GCE-REMC simulations and those from the literature.
.
Figure 2(b) shows the comparisons of binodals, the critical point of SW particles computed from isobaric GCE-REMC simulations with those from the literature more clearly. Binodals in the temperature-density plane at various pressures as listed in Table 2, displaying little deviation from those from the literature. This good agreement suggests that the novel isobaric GCE-REMC method is able to reproduce the vapor–liquid phase diagram of SW particles.
Table 2.
Table 2.
Table 2.
Vapor–liquid coexistence data of SW particles with λ = 1.5 from GCE-REMC simulations in this work. and represent the density and enthalpy of vapor binodal, and and represent the density and enthalpy of liquid binodal, respectively.
.
t
P*
0.967
0.025
131.23
0.665
–1570.53
0.02
1.067
0.053
23.80
0.612
–1410.40
0.04
1.134
0.089
–99.80
0.569
–1290.78
0.06
1.185
0.151
–261.75
0.519
–1139.45
0.08
Table 2.
Vapor–liquid coexistence data of SW particles with λ = 1.5 from GCE-REMC simulations in this work. and represent the density and enthalpy of vapor binodal, and and represent the density and enthalpy of liquid binodal, respectively.
.
Figure 3(a) shows the sampling efficiencies of various states of SW particles quantified by the enthalpy distributions at pressure close to the critical pressure . Each enthalpy distribution of 72 replicas displays the single peak, corresponding to the state in the temperature–enthalpy plane, and overlaps adjacent distributions. Even in the unstable phase-coexisting regions from H* in a range from −1000 to −500, enthalpies show unimodal distributions, rather than the multiple-peak distributions obtained from conventional canonical simulations.[47] As shown in the inset more clearly and expected according to Eq. (5), all probabilities of enthalpy follow Gaussian-like distributions and overlaps between adjacent enthalpy distributions are significant, indicating that the GCE-REMC method transforms enthalpy distributions in the phase-coexisting regions into the unimodal Gaussian-like distributions, and all possible states, including the stable or metastable vapor and liquid states, even the unstable vapor–liquid coexisting states during the vapor–liquid transition, are visited with the high efficiency.
Fig. 3. (color online) Sampling efficiencies of various states of SW particles with λ = 1.5 using the isobaric GCE-REMC method. (a) Enthalpy distributions (the left axis) of 72 replicas corresponding to the states of 72 replicas in the temperature-enthalpy plane (the right axis) at the pressure close to the critical pressure . The inset gives the details of enthalpy distribution profiles over the enthalpy ranging from −600 to −500. (b) The enthalpy distributions and temperature-enthalpy plane at pressure far from the critical pressure.
Figure 3(b) shows the enthalpy distributions of various states and temperature-enthalpy plane at the pressure far from the critical pressure. The unimodal Gaussian-like enthalpy distributions with significant overlaps between neighbor distributions also exist for stable or metastable phases, even the phase-coexisting regions, suggesting that the sampling efficiency of the GCE-REMC method is almost not affected by the variation of pressure and it still works under low pressure situations.
Figures 4(a)–4(c) present the swap acceptance rate (SAR) of 71 replica pairs using the isobaric GCE-REMC method. In this novel method, a replica exchange scheme is used to overcome the energy barrier separating different configuration regions by exchanging configurations in replica pairs with different temperatures. The exchange efficiency is characterized by the value of SAR, which measures the acceptance probability of the configuration exchanges. As shown in Figs. 4(a) and 4(b), the values of SAR at pressures and 0.04, both stay around 0.18, close to the suitable value of 0.2,[48] even for replica pairs in phase-coexisting regions. As the value of SAR depends on the overlap of enthalpy distributions in replica pairs,[49] this significant value is consistent with the evident overlap between neighbor enthalpy distributions and indicates that the configuration exchange is fulfilled with high efficiency even in phase-coexisting regions. The slight decrease for the replicas in this region is attributed to the fact that the bigger sampling slope α = 0.006 ensures that the unstable phase-coexisting region is accessed in the form of a unimodal Gaussian-like distribution, but reduces the overlaps between neighbor enthalpy distributions subsequently reducing the value of SAR. Figure 4(c) further shows the dependence of average SAR on pressure. The average value of SAR is almost a constant of 0.18 over the pressure range from 0.02 to 0.1, indicating that the efficiency of the replica exchange scheme is independent of the pressure and the GCE-REMC method works well over a wide pressure range.
Fig. 4. (color online) Swap acceptance rate (SAR) using the isobaric GCE-REMC method. ((a) and (b)) SAR of 71 replica pairs at pressures and 0.04, close to and far from the critical pressure , respectively. (c) Dependence of average value SAR of 71 replica pairs on pressure in a pressure range from 0.02 to 0.1.
4. Conclusions
In summary, we combine the isobaric replica exchange Monte Carlo (REMC) with the generalized canonical ensemble (GCE) to form a novel isobaric GCE-REMC method, and apply this method to the study of vapor–liquid transition of square-well (SW) particles. By comparing the estimated binodals and critical point with those from the literature, we find that this isobaric GCE-REMC method can reproduce the phase diagram of vapor–liquid transition. Moreover, this method shows the notable advantage that the unstable phase-coexisting regions in the vapor–liquid transition, which cannot be visited in conventional numerical simulations, are accessed with a high sampling efficiency. This new method can also work well over a wide pressure range, suggesting that it is a sophisticated tool to detect structural and thermodynamical information of metastable or unstable states during nucleation or the crystallization process in physical or biological systems.
